function [res,resEuler,p2_cup,logP2_cu] = ResidualEquationProjBondPrices(thetaVec,setup)
 
%% Getting some information from setup 
scaleX_in_g = setup.scaleX_in_g;
out         = setup.out;
x_cu        = setup.xGrid;
numX        = size(x_cu,1);
nNodes      = setup.n_nodes;
wNodes      = transpose(setup.weight_nodes);
epsNodes    = setup.epsi_nodes;
params      = setup.params;
appOrder    = setup.appOrder;
nx          = setup.nx;
ny          = setup.ny;
nx1         = setup.nx1;
gSS         = setup.gSS;
hSS         = setup.hSS;
% For the new bond price - the unknowns
thetaProj   = reshape(thetaVec,1,setup.nDim2Theta);
g0          = thetaProj(1,1); 
gx          = thetaProj(1,2:1+nx); 
if appOrder > 1 
    gxxV    = thetaProj(1,nx+2:1+nx+nx*(nx+1)/2); 
end 
if appOrder > 2 
    gxxxV   = thetaProj(1,1+nx+nx*(nx+1)/2+1:1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6); 
end 
 
%Unfold params 
BETTA= params.BETTA;
B= params.B;
CHI= params.CHI;
CHI0= params.CHI0;
THETA= params.THETA;
DELTA= params.DELTA;
ALFA= params.ALFA;
PHI= params.PHI;
PHIzero= params.PHIzero;
ZI= params.ZI;
KAPAw= params.KAPAw;
ETA= params.ETA;
RHOR= params.RHOR;
PHIpai= params.PHIpai;
PHIpai_1= params.PHIpai_1;
PHIy= params.PHIy;
PHIy_1= params.PHIy_1;
PHIc= params.PHIc;
PHIc_1= params.PHIc_1;
PHIl= params.PHIl;
NU= params.NU;
U0= params.U0;
U0d= params.U0d;
OMEGAz= params.OMEGAz;
OMEGAd= params.OMEGAd;
OMEGAn= params.OMEGAn;
OMEGAp= params.OMEGAp;
OMEGAa= params.OMEGAa;
RHOz= params.RHOz;
RHOd= params.RHOd;
RHOn= params.RHOn;
RHOp= params.RHOp;
RHOa= params.RHOa;
Kss= params.Kss;
OUTPUTss= params.OUTPUTss;
PAIss= params.PAIss;
MUZss= params.MUZss;
Css= params.Css;
AA= params.AA;
lss= params.lss;
Rss= params.Rss;
Wss= params.Wss;
 
%% Unfold the states and controls for the core of the model - all in levels 
% Getting the states at time t 
c_ba1 = out.c_ba1;
muz_cu = out.muz_cu;
d_cu = out.d_cu;
n_cu = out.n_cu;
paistar_cu = out.paistar_cu;
a_cu = out.a_cu;
 
% Getting the controls at time t 
c_cu = out.c_cu;
pai_cu = out.pai_cu;
evf_cu = out.evf_cu;
 
% Getting the states at time t+1 
c_ba1p = out.c_ba1p;
muz_cup = out.muz_cup;
d_cup = out.d_cup;
n_cup = out.n_cup;
paistar_cup = out.paistar_cup;
a_cup = out.a_cup;
 
% Getting the controls at time t+1 
c_cup = out.c_cup;
pai_cup = out.pai_cup;
evf_cup = out.evf_cup;
 
% Bond price 
p1_cup = out.p1_cup;
 
% Computing the new bond price 
x1_cu  = x_cu(:,1)/scaleX_in_g(1,1);
x2_cu  = x_cu(:,2)/scaleX_in_g(2,1);
x3_cu  = x_cu(:,3)/scaleX_in_g(3,1);
x4_cu  = x_cu(:,4)/scaleX_in_g(4,1);
x5_cu  = x_cu(:,5)/scaleX_in_g(5,1);
x6_cu  = x_cu(:,6)/scaleX_in_g(6,1);
 
p2_cu = g0(1,1)+gx(1,1).*x1_cu+gx(1,2).*x2_cu+gx(1,3).*x3_cu+gx(1,4).*x4_cu+gx(1,5).*x5_cu+gx(1,6).*x6_cu;
if appOrder > 1
    p2_cu = p2_cu+gxxV(1,1).*x1_cu.*x1_cu+2.*gxxV(1,2).*x1_cu.*x2_cu+2.*gxxV(1,3).*x1_cu.*x3_cu+2.*gxxV(1,4).*x1_cu.*x4_cu+2.*gxxV(1,5).*x1_cu.*x5_cu+2.*gxxV(1,6).*x1_cu.*x6_cu+gxxV(1,7).*x2_cu.*x2_cu+2.*gxxV(1,8).*x2_cu.*x3_cu+2.*gxxV(1,9).*x2_cu.*x4_cu+2.*gxxV(1,10).*x2_cu.*x5_cu+2.*gxxV(1,11).*x2_cu.*x6_cu+gxxV(1,12).*x3_cu.*x3_cu+2.*gxxV(1,13).*x3_cu.*x4_cu+2.*gxxV(1,14).*x3_cu.*x5_cu+2.*gxxV(1,15).*x3_cu.*x6_cu+gxxV(1,16).*x4_cu.*x4_cu+2.*gxxV(1,17).*x4_cu.*x5_cu+2.*gxxV(1,18).*x4_cu.*x6_cu+gxxV(1,19).*x5_cu.*x5_cu+2.*gxxV(1,20).*x5_cu.*x6_cu+gxxV(1,21).*x6_cu.*x6_cu;
end
if appOrder > 2
    p2_cu = p2_cu+gxxxV(1,1).*x1_cu.*x1_cu.*x1_cu+3.*gxxxV(1,2).*x1_cu.*x1_cu.*x2_cu+3.*gxxxV(1,3).*x1_cu.*x1_cu.*x3_cu+3.*gxxxV(1,4).*x1_cu.*x1_cu.*x4_cu+3.*gxxxV(1,5).*x1_cu.*x1_cu.*x5_cu+3.*gxxxV(1,6).*x1_cu.*x1_cu.*x6_cu+3.*gxxxV(1,7).*x1_cu.*x2_cu.*x2_cu+6.*gxxxV(1,8).*x1_cu.*x2_cu.*x3_cu+6.*gxxxV(1,9).*x1_cu.*x2_cu.*x4_cu+6.*gxxxV(1,10).*x1_cu.*x2_cu.*x5_cu+6.*gxxxV(1,11).*x1_cu.*x2_cu.*x6_cu+3.*gxxxV(1,12).*x1_cu.*x3_cu.*x3_cu+6.*gxxxV(1,13).*x1_cu.*x3_cu.*x4_cu+6.*gxxxV(1,14).*x1_cu.*x3_cu.*x5_cu+6.*gxxxV(1,15).*x1_cu.*x3_cu.*x6_cu+3.*gxxxV(1,16).*x1_cu.*x4_cu.*x4_cu+6.*gxxxV(1,17).*x1_cu.*x4_cu.*x5_cu+6.*gxxxV(1,18).*x1_cu.*x4_cu.*x6_cu+3.*gxxxV(1,19).*x1_cu.*x5_cu.*x5_cu+6.*gxxxV(1,20).*x1_cu.*x5_cu.*x6_cu+3.*gxxxV(1,21).*x1_cu.*x6_cu.*x6_cu+gxxxV(1,22).*x2_cu.*x2_cu.*x2_cu+3.*gxxxV(1,23).*x2_cu.*x2_cu.*x3_cu+3.*gxxxV(1,24).*x2_cu.*x2_cu.*x4_cu+3.*gxxxV(1,25).*x2_cu.*x2_cu.*x5_cu+3.*gxxxV(1,26).*x2_cu.*x2_cu.*x6_cu+3.*gxxxV(1,27).*x2_cu.*x3_cu.*x3_cu+6.*gxxxV(1,28).*x2_cu.*x3_cu.*x4_cu+6.*gxxxV(1,29).*x2_cu.*x3_cu.*x5_cu+6.*gxxxV(1,30).*x2_cu.*x3_cu.*x6_cu+3.*gxxxV(1,31).*x2_cu.*x4_cu.*x4_cu+6.*gxxxV(1,32).*x2_cu.*x4_cu.*x5_cu+6.*gxxxV(1,33).*x2_cu.*x4_cu.*x6_cu+3.*gxxxV(1,34).*x2_cu.*x5_cu.*x5_cu+6.*gxxxV(1,35).*x2_cu.*x5_cu.*x6_cu+3.*gxxxV(1,36).*x2_cu.*x6_cu.*x6_cu+gxxxV(1,37).*x3_cu.*x3_cu.*x3_cu+3.*gxxxV(1,38).*x3_cu.*x3_cu.*x4_cu+3.*gxxxV(1,39).*x3_cu.*x3_cu.*x5_cu+3.*gxxxV(1,40).*x3_cu.*x3_cu.*x6_cu+3.*gxxxV(1,41).*x3_cu.*x4_cu.*x4_cu+6.*gxxxV(1,42).*x3_cu.*x4_cu.*x5_cu+6.*gxxxV(1,43).*x3_cu.*x4_cu.*x6_cu+3.*gxxxV(1,44).*x3_cu.*x5_cu.*x5_cu+6.*gxxxV(1,45).*x3_cu.*x5_cu.*x6_cu+3.*gxxxV(1,46).*x3_cu.*x6_cu.*x6_cu+gxxxV(1,47).*x4_cu.*x4_cu.*x4_cu+3.*gxxxV(1,48).*x4_cu.*x4_cu.*x5_cu+3.*gxxxV(1,49).*x4_cu.*x4_cu.*x6_cu+3.*gxxxV(1,50).*x4_cu.*x5_cu.*x5_cu+6.*gxxxV(1,51).*x4_cu.*x5_cu.*x6_cu+3.*gxxxV(1,52).*x4_cu.*x6_cu.*x6_cu+gxxxV(1,53).*x5_cu.*x5_cu.*x5_cu+3.*gxxxV(1,54).*x5_cu.*x5_cu.*x6_cu+3.*gxxxV(1,55).*x5_cu.*x6_cu.*x6_cu+gxxxV(1,56).*x6_cu.*x6_cu.*x6_cu;
end
 
% The bond price next period 
if setup.computeP2_cup == 1 
    xVec_cup = zeros(numX,nx);
    mx = setup.mx;
    if mx > 0 
    end 
    if setup.myx > 0 
        xVec_cup(:,mx+1) = c_cu(:,1)-gSS(1,1);
    end 
    diaghx = diag(setup.hx);
    xVec_cup(:,nx1+1:nx)  = x_cu(:,nx1+1:nx).*repmat(transpose(diaghx(nx1+1:nx)),numX,1); 
    
    %Adding shocks to the states at t + 1
    diagEta = diag(setup.eta(nx1+1:nx,:));
    x1_cup  = repmat(xVec_cup(:,1),1,nNodes);
    x2_cup  = repmat(xVec_cup(:,2),1,nNodes) +2./(exp(-OMEGAz.*x_cu(:,2))+1).*diagEta(1,1).*repmat(epsNodes(1,:),numX,1);
    x3_cup  = repmat(xVec_cup(:,3),1,nNodes) +2./(exp(-OMEGAd.*x_cu(:,3))+1).*diagEta(2,1).*repmat(epsNodes(2,:),numX,1);
    x4_cup  = repmat(xVec_cup(:,4),1,nNodes) +2./(exp(-OMEGAn.*x_cu(:,4))+1).*diagEta(3,1).*repmat(epsNodes(3,:),numX,1);
    x5_cup  = repmat(xVec_cup(:,5),1,nNodes) +2./(exp(-OMEGAp.*x_cu(:,5))+1).*diagEta(4,1).*repmat(epsNodes(4,:),numX,1);
    x6_cup  = repmat(xVec_cup(:,6),1,nNodes) +2./(exp(-OMEGAa.*x_cu(:,6))+1).*diagEta(5,1).*repmat(epsNodes(5,:),numX,1);
 
    % Scaling the states at time t+1 by the stdX
    x1_cup  = x1_cup/scaleX_in_g(1,1);
    x2_cup  = x2_cup/scaleX_in_g(2,1);
    x3_cup  = x3_cup/scaleX_in_g(3,1);
    x4_cup  = x4_cup/scaleX_in_g(4,1);
    x5_cup  = x5_cup/scaleX_in_g(5,1);
    x6_cup  = x6_cup/scaleX_in_g(6,1);
 
    p2_cup = g0(1,1)+gx(1,1).*x1_cup+gx(1,2).*x2_cup+gx(1,3).*x3_cup+gx(1,4).*x4_cup+gx(1,5).*x5_cup+gx(1,6).*x6_cup;
    if appOrder > 1
        p2_cup = p2_cup+gxxV(1,1).*x1_cup.*x1_cup+2.*gxxV(1,2).*x1_cup.*x2_cup+2.*gxxV(1,3).*x1_cup.*x3_cup+2.*gxxV(1,4).*x1_cup.*x4_cup+2.*gxxV(1,5).*x1_cup.*x5_cup+2.*gxxV(1,6).*x1_cup.*x6_cup+gxxV(1,7).*x2_cup.*x2_cup+2.*gxxV(1,8).*x2_cup.*x3_cup+2.*gxxV(1,9).*x2_cup.*x4_cup+2.*gxxV(1,10).*x2_cup.*x5_cup+2.*gxxV(1,11).*x2_cup.*x6_cup+gxxV(1,12).*x3_cup.*x3_cup+2.*gxxV(1,13).*x3_cup.*x4_cup+2.*gxxV(1,14).*x3_cup.*x5_cup+2.*gxxV(1,15).*x3_cup.*x6_cup+gxxV(1,16).*x4_cup.*x4_cup+2.*gxxV(1,17).*x4_cup.*x5_cup+2.*gxxV(1,18).*x4_cup.*x6_cup+gxxV(1,19).*x5_cup.*x5_cup+2.*gxxV(1,20).*x5_cup.*x6_cup+gxxV(1,21).*x6_cup.*x6_cup;
    end
    if appOrder > 2
        p2_cup = p2_cup+gxxxV(1,1).*x1_cup.*x1_cup.*x1_cup+3.*gxxxV(1,2).*x1_cup.*x1_cup.*x2_cup+3.*gxxxV(1,3).*x1_cup.*x1_cup.*x3_cup+3.*gxxxV(1,4).*x1_cup.*x1_cup.*x4_cup+3.*gxxxV(1,5).*x1_cup.*x1_cup.*x5_cup+3.*gxxxV(1,6).*x1_cup.*x1_cup.*x6_cup+3.*gxxxV(1,7).*x1_cup.*x2_cup.*x2_cup+6.*gxxxV(1,8).*x1_cup.*x2_cup.*x3_cup+6.*gxxxV(1,9).*x1_cup.*x2_cup.*x4_cup+6.*gxxxV(1,10).*x1_cup.*x2_cup.*x5_cup+6.*gxxxV(1,11).*x1_cup.*x2_cup.*x6_cup+3.*gxxxV(1,12).*x1_cup.*x3_cup.*x3_cup+6.*gxxxV(1,13).*x1_cup.*x3_cup.*x4_cup+6.*gxxxV(1,14).*x1_cup.*x3_cup.*x5_cup+6.*gxxxV(1,15).*x1_cup.*x3_cup.*x6_cup+3.*gxxxV(1,16).*x1_cup.*x4_cup.*x4_cup+6.*gxxxV(1,17).*x1_cup.*x4_cup.*x5_cup+6.*gxxxV(1,18).*x1_cup.*x4_cup.*x6_cup+3.*gxxxV(1,19).*x1_cup.*x5_cup.*x5_cup+6.*gxxxV(1,20).*x1_cup.*x5_cup.*x6_cup+3.*gxxxV(1,21).*x1_cup.*x6_cup.*x6_cup+gxxxV(1,22).*x2_cup.*x2_cup.*x2_cup+3.*gxxxV(1,23).*x2_cup.*x2_cup.*x3_cup+3.*gxxxV(1,24).*x2_cup.*x2_cup.*x4_cup+3.*gxxxV(1,25).*x2_cup.*x2_cup.*x5_cup+3.*gxxxV(1,26).*x2_cup.*x2_cup.*x6_cup+3.*gxxxV(1,27).*x2_cup.*x3_cup.*x3_cup+6.*gxxxV(1,28).*x2_cup.*x3_cup.*x4_cup+6.*gxxxV(1,29).*x2_cup.*x3_cup.*x5_cup+6.*gxxxV(1,30).*x2_cup.*x3_cup.*x6_cup+3.*gxxxV(1,31).*x2_cup.*x4_cup.*x4_cup+6.*gxxxV(1,32).*x2_cup.*x4_cup.*x5_cup+6.*gxxxV(1,33).*x2_cup.*x4_cup.*x6_cup+3.*gxxxV(1,34).*x2_cup.*x5_cup.*x5_cup+6.*gxxxV(1,35).*x2_cup.*x5_cup.*x6_cup+3.*gxxxV(1,36).*x2_cup.*x6_cup.*x6_cup+gxxxV(1,37).*x3_cup.*x3_cup.*x3_cup+3.*gxxxV(1,38).*x3_cup.*x3_cup.*x4_cup+3.*gxxxV(1,39).*x3_cup.*x3_cup.*x5_cup+3.*gxxxV(1,40).*x3_cup.*x3_cup.*x6_cup+3.*gxxxV(1,41).*x3_cup.*x4_cup.*x4_cup+6.*gxxxV(1,42).*x3_cup.*x4_cup.*x5_cup+6.*gxxxV(1,43).*x3_cup.*x4_cup.*x6_cup+3.*gxxxV(1,44).*x3_cup.*x5_cup.*x5_cup+6.*gxxxV(1,45).*x3_cup.*x5_cup.*x6_cup+3.*gxxxV(1,46).*x3_cup.*x6_cup.*x6_cup+gxxxV(1,47).*x4_cup.*x4_cup.*x4_cup+3.*gxxxV(1,48).*x4_cup.*x4_cup.*x5_cup+3.*gxxxV(1,49).*x4_cup.*x4_cup.*x6_cup+3.*gxxxV(1,50).*x4_cup.*x5_cup.*x5_cup+6.*gxxxV(1,51).*x4_cup.*x5_cup.*x6_cup+3.*gxxxV(1,52).*x4_cup.*x6_cup.*x6_cup+gxxxV(1,53).*x5_cup.*x5_cup.*x5_cup+3.*gxxxV(1,54).*x5_cup.*x5_cup.*x6_cup+3.*gxxxV(1,55).*x5_cup.*x6_cup.*x6_cup+gxxxV(1,56).*x6_cup.*x6_cup.*x6_cup;
    end
else
    p2_cup = NaN;
end 
 
% The Euler errors 
E_sdf_p1_cup = (BETTA.*exp(-d_cu).*exp(-pai_cup).*exp(d_cup).*exp(muz_cup).^(CHI.*(CHI0 - 1) - CHI0).*(exp(c_cu) - B.*exp(-muz_cu).*exp(c_ba1)).^CHI.*(AA./(exp(evf_cu).^(1./(ALFA - 1)).*exp(muz_cup).^((CHI - 1).*(CHI0 - 1)).*(exp(d_cup).*(((exp(c_cup) - B.*exp(-muz_cup).*exp(c_cu))./Css.^CHI0).^(1 - CHI)./(CHI - 1) - U0d + (PHIzero.*exp(n_cup).*(1 - 1./(-(exp(-a_cup).*(exp(c_cup) + DELTA.*Kss))./(Kss.^THETA.*((ZI.*(exp(pai_cup)./PAIss.^NU - 1).^2)./2 - 1))).^(1./(THETA - 1))).^(1 - 1./PHI))./(1./PHI - 1)) - U0 + (AA.*BETTA)./exp(evf_cup).^(1./(ALFA - 1))))).^ALFA)./(exp(c_cup) - B.*exp(-muz_cup).*exp(c_cu)).^CHI.*exp(p1_cup);
f1 = E_sdf_p1_cup./exp(p2_cu) - 1;
% Numerical integration 
resEuler = f1*wNodes;
logP2_cu = log(E_sdf_p1_cup*wNodes);
resEulerWeight = resEuler.*(1./(1+sum(x_cu.^2,2)));
res = resEulerWeight(:);
end
